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Abstract 

Unsteady Reynolds-averaged Navier-Stokes (RANS) 
computations are presented for the flow past a flapping 
NACA 0012 aerofoil to study the effect of various motion 
parameters on the propulsive efficiency and thrust 
coefficient. The RANS solver used to obtain the time- 
accurate solutions is based on an implicit dual time stepping 
scheme with finite volume nodal point spatial discretization. 
The algebraic eddy viscosity model due to Baldwin and 
Lomax is used for the turbulence closure. The motion 
parameters considered are reduced frequency, amplitude of 
the pitch and plunge motions and phase shift between them. 
The computed results compare well with the available data 
in the literature. 
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Introduction 

In recent years, Micro Air Vehicles (MAVs), 
Unmanned Air Vehicles (UAV's) and Nano Air 
Vehicles (NAV's) are becoming increasingly important, 
especially in the area of military/defence surveillance. 
These are classified into fixed, rotary and flapping 
wing air vehicles. At low range Reynolds number, 
flapping wings are more efficient and more easily 
manoeuvrable compared to fixed wings. Numerous 
computational as well as experimental studies have 
been conducted to investigate the kinematics and 
dynamics of flapping wings. 

Many important features of flapping aerofoil 
behaviour are depicted by the classical linear theory. 
Theodorsen developed compact expressions for forces 
and moments of a flapping aerofoil under the 
assumption of small perturbation for inviscid and 
incompressible flow. The flow is divided into two 
classes: the non-circulating flow and the circulatory 



flow due to wake vortices. The thrust force and the 
power input experienced by the flapping aerofoil were 
given by Garrick . Ho and Chen studied the unsteady 
wake of a plunging Aerofoil NACA 0012 in a low 
turbulence wind- tunnel by a miniature multiple hot- 
wire probe at Reynolds number (2.1 * 10 4 ~ 10 5 ), 
reduced frequency (0 to 1), and mean angle of attack 5°. 
Under the condition of no dynamic stall, the near 
wake consists of viscous and inviscid parts. The 
viscous wake is, the result of the merging of two 
boundary layers on the aerofoil, has high velocity 
gradient and turbulence levels, and is limited to a very 
thin region. The inviscid wake is caused by the 
induced flow of the circulation around the plunging 
aerofoil and is laminar and of large width. Anderson 
and Anderson et al. have obtained visualization and 
force data for a plunging and pitching aerofoil moving 
in a water tank facility over a large range of 
amplitudes, frequencies, and phase angles. Propulsive 
efficiency, as high as 87% is measured experimentally 
under conditions of optimal wake formation. They 
found that agreement between the experimental data 
and numerical predictions of a non linear 
incompressible unsteady potential-flow method is 
good when either very weak or no flow separation 
vortices forms. Isogai et al. simulated dynamic stall 
phenomena around a flapping NACA 0012 aerofoil by 
a Navier-Stokes code at a free stream Mach number of 
0.3 and Reynolds number of 10 5 . The Baldwin and 
Lomax algebraic turbulence model is used in the 
computation. They found that high propulsive 
efficiency occurs for the case in which the pitching 
oscillation advances 90 degree ahead of the plunging 
oscillation and the reduced frequency is at some 
optimum value, for which there appears no 
appreciable flow separation in spite of large-amplitude 
oscillations. 
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IMPRANS Solver 

The solver is based on an implicit finite volume nodal 
point spatial discretization scheme with dual time 
stepping. Inviscid flux vectors are calculated by using 
the flow variables at the six neighbouring points of 
hexahedral volume. Turbulence closure is achieved 
through the algebraic eddy viscosity model of Baldwin 
and Lomax . 

The two-dimensional Reynolds averaged Navier- 
Stokes equations for a moving domain can be written 
in non-dimensional conservative form as 



dU dF dG 

dt dx dy 



dV dW 

dx dy 



(1) 

Here x and y are the Cartesian coordinates and f is the 

time while x and ^ are the Cartesian velocity 
components of the moving domain. For a fixed 

domain, the gird speeds x and ^ are zero. U is the 
vector of conserved variables; F, G are inviscid flux 
vectors and V, W are viscous flux vectors. 

Applying Euler's implicit time differencing formula 
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to the governing Eq. (1), we obtain 
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Here U" = U (t)= U (n At) is the solution vector at time 
level n and AU" = (U n+1 - U") is the change in U"over 
time step At. In order to facilitate the finite volume 
formulation, the above equations are written in the 
integral form as 



\\ a AU n dxdy + At\ T ( F _ V f +l dy _ ( G _ W ) n+l dx 



(4) 



where O is any two-dimensional flow domain and f 
is the boundary curve. 

In the nodal point finite volume approach, the flow 
variables are associated with each mesh point of the 
grid and the integral conservative equations are 
applied to each control volume obtained by joining the 
centroids of the four neighbouring cells of a nodal 
point. 



Application of nodal point spatial discretization to 
Eq. (4) leads to the following equations for the 
computational cell Oij 



U L 
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= 



(5) 



where hj is the area of quadrilateral and the integral 
refers to a contour integration around the boundary Fij 
of the cell Qij in anticlockwise direction. The fluxes are 
calculated across the four sides of the control volume. 

Linearising the changes in flux vectors using Taylor's 
series expansions in time and assuming locally 
constant transport properties, Eq. (5) can be simplified 
to 
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Here A, B, R and S are the Jacobian matrices given by 
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The algebraic eddy viscosity model due to Baldwin 
and Lomax is used for turbulence closure. This RANS 
solver has been extensively validated for computing 
unsteady flow past pitching aerofoils and wings, 
plunging aerofoils, pitching and plunging aerofoils 
helicopter rotor blades, wind turbines etc. Here, the 
solver has been applied to study the effect of various 
parameters on flapping motion of NACA 0012 aerofoil 
and compared with available results of Isogai et al. 
and Tuncer et al. 

Results 

Here, NACA 0012 aerofoil is undergoing combined 
pitching and plunging motion sinusoidally is 
considered at different motion parameters. The C-H 
grid of size 247 * 65 around the aerofoil has been 
generated using commercial software GRIDGEN. The 
grid points are properly clustered near the leading 
edge, trailing edge and along wall normal direction 
which is shown in Fig. 1(a) and (b). 
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FIG. 1(a) C-H GRID AROUND NACA0012 AEROFOIL 




FIG. 1(b) CLOSE UP VIEW OF GRID 

The combined pitching and plunging oscillation of an 
aerofoil is shown in Fig. 2. The plunging motion is 
defined as the vertical motion at right angles to the 
direction of uniform free stream velocity The 
NACA 0012 aerofoil with chord length c performs a 
sinusoidal plunging. The position of the aerofoil is y (f) 
given by 

y{t) = yosin{wt) (8) 
where f is physical time, oj the angular frequency. 




FIG. 2 AEROFOIL IN COMBINED PITCHING AND PLUNGING 
MOTION WITH A PHASE DIFFERENCE OF = 90° 

Plunge amplitude, yo, is positive in the upward 
direction and the reduced frequency is given by, 

k = aclU 



The instantaneous non-dimensional plunging velocity 
is 



ylU„ = h a kcos(oat) 



(9) 



where the dot denotes a differentiation with respect to 
t and the non-dimensional plunge amplitude is 

ha = y Ic. 

The instantaneous effective angle of attack due to pure 
plunging is 



a = a sm(a>t + < 



(10) 



is 



where the amplitude of pitching oscillation is a Q , 
the phase angle ahead of the plunging motion which is 
also shown in Fig. 2. The instantaneous lift and thrust 
coefficients are C; and Ct, respectively. The 
instantaneous pitching moment coefficient around the 
pitching-pivot point is Cm (positive in the nose-up 
sense). 



The instantaneous input power co-efficient is 



C„ 



-ac,y + C m cd)/U,, 



(11) 



The time-averaged lift coefficient, pitching moment 
coefficient, thrust coefficient and input power 
coefficient over an oscillation period are G, Cm, Ct and 
Cp respectively. The mean thrust coefficient is 
therefore defined as 



(c,)=-c„ +(C d )„ 



(12) 



where Ca is the mean drag coefficient, averaged for one 
flapping period. (Q)steadyis the steady drag of the non- 
moving wing at its present mean angle of attack. 

The propulsive efficiency is calculated from the ratio 
between power output and power input. 
If dimensionless coefficients are used, this is equal to 
the ratio of mean thrust coefficient to mean power 
input coefficient 

kJ=(c,)/(cJ (13) 

As it can be identified from Fig. 2, the combined 
motion produces wave motion propagating in the x- 
direction. It is clear that the origin of this wave motion 
can be attributed to the coupled heaving and pitching 
oscillations with a finite phase difference. When we 
define the propagation velocity of the wave, the ratio 
of wave velocity to free stream velocity depends on 
the critical reduced frequency which in turn is a 
function of pivot point a, aol ha and cp. However, the 
behaviour of critical reduced frequency is still useful 
to select the parameters such as k, cp, ha and ao for 
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which the numerical simulations are to be performed 
for various combinations of k and cp. For this reason, 
the behaviour of critical reduced frequency is studied 
or simulated for the following two different plunge 
amplitudes (ha), 

1. ha= 1.0, ao=20°, a = 0, M = 0.3 and Re = 1.0 * 10 5 

2. h>=2.0, ao= 10°, a = 0, M = 0.3 and Re= 1.0 * 10 s 

The present numerical simulations have been 
performed for various combinations of k and & for 
each of these two cases. 

1. ha = 1.0, ao=20°, a = 0,M = 0.3 and Re = 1.0 x 10 s 

Table 1 shows the results optimized for maximum 
propulsion efficiency, r/p and time-averaged thrust 
coefficient, C t, respectively with respect to k for 
cp = 90°, which have been compared with Navier- 
Stokes results available in Isogai et al. The highest 
propulsion efficiency of 79.8% with a time-averaged 
thrust coefficient of 0.310 and the highest time- 
averaged thrust coefficient of 0.6118 with propulsion 
efficiency 43.10% are obtained. Fig. 3 shows the 
coefficient of lift, drag and moment versus the angle of 
attack for the final converged cycle. The pressure 
contours are shown in Fig. 4. 

TABLE 1 COEFFICIENT OF THRUST AND PROPULSION 
EFFICIENCY FOR A PITCHING-PLUNGING NACA 0012 
AEROFOIL 



k 


Coefficient of thrust 


Propulsive efficiency x 
100 


Present 
(RANS) 


Isogai [6] 


Present 
(RANS) 


Isogai [6] 


0.5 


0.10090528 


0.310 


0.79840 


0.725 


0.6 


0.20012392 


0.466 


0.78002 


0.6994 


0.7 


0.31576321 


0.524 


0.75058 


0.6742 


0.8 


0.44585174 


0.566 


0.70883 


0.6205 


0.9 


0.51040155 


0.595 


0.53349 


0.5815 


1.0 


0.57188636 


0.6118 


0.43102 


0.500 



2. ha = 2.0,ao=10°,a = 0,M = 0.3 and Re = 1.0 * 10 s 

Table 2 and Table 3 show propulsion efficiency (t]p) 
and time-averaged thrust coefficient (C t) respectively 
which are compared against cp with k as a varying 
parameter compared to the available results obtained 
by Isogai et al. and Tuncer et al. both of who used a 
Navier-Stokes code with Baldwin and Lomax 
turbulence model. Although the agreements between 



the two computations are fairly good for cp > 60°, the 
large quantitative discrepancies are seen in both r/p 
and Cf for cp = 30°, in which case, the severe leading 
edge-flow separations have been observed in our 
computations for all of the reduced frequencies 
computed, whereas no flow separation has been 
observed for k = 0.15 in the computations by Tuncer 
and Platzer. It is interesting to see their results that the 
highest C t is obtained for k = 0.15 and cp = 60°. It is 
clear that the low values of C t and poor efficiency 
obtained for k = 0.5 can be attributed to the occurrence 
of large-scale leading-edge separation. The highest 
propulsion efficiency of 95.14% with a time-averaged 
thrust coefficient of 0.0875 and the highest time- 
averaged thrust coefficient of 0.3168 with propulsion 
efficiency 19.09% are obtained. Fig. 5 shows the 
coefficient of lift, drag and moment versus the angle of 
attack for the final converged cycle. The Mach number 
contour at different instants of time for one complete 
cycle of flapping motion of the aerofoil is plotted in 
Figs. 6. The critical flow patterns are well captured and 
flow separations are seen clearly in the plots. There is 
a difference in predicted propulsion efficiency and 
time-averaged thrust coefficient values and that of 
Isogai et al. at smaller reduced frequencies is probably 
due to the presence of viscous effects. 

TABLE 2 COEFFICIENT OF THRUST FOR A PITCHING- 
PLUNGING NACA 0012 AEROFOIL 



V 


k 


Coefficient of thrust 


Present 
RANS 


Isogai [6] 


Tuncer [16] 


30° 


0.15 


0.1267 


0.1 


1.4 


0.3 


0.01601 


0.03 


0.32 


0.5 


0.0178 


0.04 


0.1 


60° 


0.15 


0.112166 


0.71 


1.325 


0.3 


0.131885 


0.5 


0.5989 


0.5 


0.17113 


0.26 


0.34 


90° 


0.15 


0.08750 


0.88 


1.015 


0.3 


0.21918 


0.706 


0.752 


0.5 


0.28837 


0.415 


0.41 


120° 


0.15 


0.08895 


0.66 


0.72 


0.3 


0.212377 


0.74 


0.725 


0.5 


0.316890 


0.409 


0.389 


150° 


0.15 


0.1061 


0.665 


0.66 


0.3 


0.14267 


0.4401 


0.492 


0.5 


0.25199 


0.26 


0.2599 
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TABLE 3 PROPULSION EFFICIENCY FOR A PITCHING- 
PLUNGING NACA 0012 AEROFOIL 



<p 


k 


Propulsive efficiency x 100 


RANS 




i unccr L-ioj 


30° 


0.15 


0.8098 


0.05 


0.64 


0.3 


03405 


0.036 


0.15 


0.5 




0.05 


0.089 


60° 


0.15 


931473 


0.556 


0.72 


0.3 


0.271101 


0.2722 


0.32 


0.5 


0.11772 


0.142 


0.164 


90° 


0.15 


0.95145 


0.787 


0.832 


0.3 


0.40189 


0.406 


0.432 


0.5 


0.18597 


0.215 


0.22 


120° 


0.15 


0.81822 


0.765 


0.865 


0.3 


0.37524 


0.391 


0.42 


0.5 


0.19096 


0.198 


0.201 


150° 


0.15 


0.6952 


0.62 


0.72 


0.3 


0.2024 


0.1982 


0.2535 


0.5 


0.1441 


0.1 


0.152 




» 0.-20° 
■ 831TO y/c = -0.00027 

S2?r»+ 


a= 1SCE° 
y/c= 0.3824 


a- 12.47° 
y/c=0.7816 


tt- -4.45° 
y/c = 9749 


a.- -18.01° 
y/c = 0.434 


11.-181° 
y/c = -0.4335 


at« -4 45 ° 
y/c = -0.9748 


a - 12.464° 
y/c= -0.782 


(1-18.01° 
y/c = -0.434 


a=20° 
y/c=-0 00034 



FIG. 4 THE PRESSURE FIELDS AT DIFFERENT INSTANTS OF 
TIME FOR ONE CYCLE OF FLAPPING MOTION OF AEROFOIL 
AT fc,=1.0, ao= 20°, a = 1/2, M = 0.3, k = 0.9 and <p = 90° 




FIG. 3 THE COEFFICIENT OF LIFT, DRAG AND MOMENT 
VERSUS THE ANGLE OF ATTACK FOR THE FINAL CYCLE AT 
h,= 1.0, ao= 20°, a = 1/2, M = 0.3, k = 0.9 and cp = 90° 



FIG. 5 THE COEFFICIENT OF LIFT, DRAG AN D MOMENT 
VERSUS THE ANGLE OF ATTACK FOR THE FINAL CYCLE AT 
ha= 2.0, ao= 10°, a = 1/2, M = 0.3, k = 0.5 and <p = 30° 
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q J ices 
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a = 9.435° 
y/c= 1.3297 


0311173 
0JJ1139JS 

« 0.111+14 | 
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, . ... 


a, = 4.2879° 
■ y/c= 1.646 




* 


tx= -2.781 ° 
y/c = 0.4529 


y/c = -119292 
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a. — 9.8706° 
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a =6.325° 
y/c=-1.86 
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y/c= -0.9296 


3=5° 
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FIG. 6: THE MACH CONTOURS AT DIFFERENT INSTANTS OF 
TIME FOR ONE CYCLE OF FLAPPING MOTION OF AEROFOIL 
AT h„. 2.0, ao= 10°, a = 1/2, M = 0.3, k = 0.5 and <p = 30° 



Conclusions 

The time-averaged thrust coefficient of a combined 
pitching-plunging NACA 0012 aerofoil computed by 
the RANS solver agrees fairly well with those obtained 
by Navier-Stokes solutions in the literature. The effect 
of various motion parameters on the time-averaged 
thrust generation and propulsive efficiency has been 
studied. The highest propulsive efficiency and the 
highest thrust coefficient do not occur at the same 
reduced frequency. Higher propulsive efficiency 
usually occurs at lower reduced frequency, where as 
higher thrust occurs at higher reduced frequency. As 
the flow becomes more and more unsteady with 
increasing reduced frequency, a large amount of 
vorticity is shed from the trailing edge. At higher 
values of pitching amplitude time histories of the 
coefficient of lift, thrust, drag and moment are smooth 
purely sinusoidal but, as the pitching amplitude is 
decreased, the time histories of these coefficients are 
not sinusoidal as their shape is deformed thus 
resulting in the decrease in propulsion efficiency and 
an increase in the thrust coefficient comparatively. 
The maximum propulsive efficiency and minimum 



thrust for the range of cases considered in the present 
study occur when the phase angle between pitching 
and plunging motion is approximately 90°. 
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